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We use analytic and numerical methods to analyze the dynamics of vortices following the quench of 
a Type-II superconductor under the application of an external magnetic field. In three dimensions, 
in the absence of a field, the spacing between vortices scales with time t with an exponent (f> = 
0.414 ± 0.01, In a thin sheet of superconductor, the scaling exponent is (f> — 0.294 ± 0.01. When an 
external magnetic field h is applied, the vortices are confined with respect to the length scale of the 
Abrikosov lattice, leading to a crossover between the power-law scaling length scale and the lattice 
length scale. From this we suggest a one-parameter scaling of r with h and r that is consistent with 
numerical data. 
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I. INTRODUCTION 

Systems with continuous symmetries exhibit interest- 
ing time-dependent scaling behavior following a quench, 
due to the dynamics of topological defects 0,0,0,0,010 
7]. Physical realizations include superfluid He 4 and He 3 , 
Bose-Einstein condensates, liquid crystals, and models 
of cosmic strings - defects produced during phase transi- 
tions in the early universe 0]. In this paper, we investi- 
gate the effect of topological defects in the approach to 
equilibrium of Type-II superconductors and examine the 
consequences to the dynamics of the presence of a gauge 
field HEl. 

Following a quench from the normal phase, a Type-II 
superconductor is threaded by a series of vortex cores 
(which are line defects) due to fluctuations that had ex- 
isted in the normal phase. Over time, these vortex cores 
interact through local and long-range forces. In the ab- 
sence of an external field, the number of vortices de- 
creases with time, giving rise to a characteristic length- 
scale r(t) that coarsens with time t, generically with a 
power law whose exponent we denote by <f>. 

A variety of predictions have been made for the scal- 
ing of r(t) by treating the problem with various reduced 
models. A first approximation to the scaling of defect 
spacing can be found by solving the equations of motion 
for a pair of defects. For local interactions between de- 
fects the strongest contribution is from the immediate 
neighbor, and the effect of the remainder of the defects 
is to distribute momentum via long-wavelength modes 
throughout the bulk, leading to over-damped dissipative 
dynamics. 

We can estimate the dynamics in this approximation 
from the potential energy of interaction of a defect pair 
0. This will play a role in determining how clusters 
of defects expand, which in turn determines the time- 
dependence of such quantities as the total defect volume, 
spacing between defects, and defect correlation functions. 
If the interaction is purely from the stress induced by 
the topological defect in the continuously symmetric field 



then the potential will be of the form V(r) oc log(r/£) for 
both point defects in 2D and line defects in 3D, leading 
to a force of the form F(r) oc (r/^)" 1 . Here we have 
introduced a characteristic length scale £, which is of or- 
der the vortex core size. Similarly, a point defect in 3D 
would experience a potential V(r) oc {r/C)~ l an d a force 
proportional to . 

In overdamped dynamics, the time-derivative r(t) is 
proportional to the force of interaction between two de- 
fects. For the r _1 force appropriate for line defects in 
three dimensions, we find that when r(t) is much greater 
than the initial inter-vortex distance, r oc with the 
dynamical exponent (f> = 1/2. 

In a superconductor, vortices are surrounded by loops 
of superconducting currents, matched to a magnetic field 
which threads the vortices. As a result, there are mag- 
netic interactions which depend on the geometry and di- 
mension of the system in addition to the direct interac- 
tion of the topological defects. In the bulk of a super- 
conductor the field is screened and the force law is ex- 
ponentially decaying with the distance, leading to a log- 
arithmic scaling of the inter-vortex distance with time. 
In a thin film of superconducting material (which will 
be referred to as the 2DF system), where the magnetic 
flux is allowed to penetrate into the third dimension, the 
force law is that of point interactions in three dimensional 
space, and results in a power law scaling of r(t) with an 
exponent <p = 1/3. 

These exponents are close to those observed in experi- 
ment and simulation, but not exactly so. As first empha- 
sized by Toussaint and Wilczek 0, fluctuations in the 
initial conditions can lead to a distribution of topological 
defects that may be long-lived, depending on the spatial 
dimension, thus influencing the long time scaling expo- 
nent 4>. In the context of vortices, such effects lead 
to the prediction that = 3/8, with subsequent refine- 
ments [l(j, [ll| yielding the predictions <f> — 3/10 for the 
2DF system and <j> = 3/7 for the fully 3D case. 

The purpose of this paper is to test these predictions 
and extend them to the case of non-zero external applied 
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Citation 


System 





Ref [12] 
Ref [£] 
Ref [3] 
Ref [10] 
Ref [llj 


Complex TDLG 

3DXY model 
Complex TDLG 
Langevin dynamics 
2D Monte-Carlo 


0.375 ± 0.03 
0.44 ±0.01 
0.45 ±0.01 
0.45 ± 0.05 
0.42 ± 0.02 



TABLE I: Literature predictions for the dynamic scaling ex- 
ponent 4> 



magnetic field. Our simulations of the 2DF and 3D bulk 
superconductor systems yield sca ling laws that are con- 
sistent with the predictions of [lfj, . In the presence of 
an externally applied magnetic field h (where h is in units 
of the critical field H c ), we found that the vortex dynam- 
ics approach predicts a scaling form for r(t, h), and that 
our numerical data collapse to a universal curve after the 
initial period of vortex annihilation. 

A summary of the various predictions for the dynamic 
scaling exponent is given in Table [Q 

This paper is organized as follows. In section[H]wc give 
in more detail the argument to derive the scaling laws 
from vortex dynamics in all three systems taking into ac- 
count the effect of the fields generated by superconduct- 
ing currents. We describe in section ITTll an improvement 
on the vortex dynamic calculations which has been used 
to make a more accurate prediction of the vortex scaling 
laws 0, El El ■ This method takes into account vortex- 
vortex annihilation and gives the observed scaling laws 
for both zero and non-zero external magnetic fields. For 
the latter case, we suggest a scaling form for the vortex 
separation r as a function of h and t such that the data 
collapse onto one universal curve: 



a/XVh-r(t)/X = f(texp(-a/xVh)), 



(1) 



where a = X^/<Po, <&o is the quantized magnetic flux of 
a single vortex, and A is the penetration depth of the 
magnetic field into the superconducting bulk. In section 
IIVI we describe the construction of a simulation to test 
these scaling laws, and compare in section the results 
of simulation with the predicted scaling and magnetic 
field data collapse. We observe the scalings predicted by 
the vortex annihilation models, and the non-zero external 
field data is collapsed onto a single curve by our proposed 
choice of variables. 



II. OVERDAMPED VORTEX DYNAMICS 

We now proceed to estimate the scaling laws for the 
inter- vortex separation by considering the dynamics of a 
pair of vortices as representative of the full many-body 
system-wide behavior. There are forces between nearby 
vortices induced both by the effect they have on the phase 
of the wavefunction, and by the magnetic flux which 



threads them. The interaction through the phase field 
falls off at a different rate than the magnetic interaction, 
and so there are multiple scaling regimes. 

We earlier gave a heuristic argument for the 0=1/2 
scaling law arising from the phase field interactions. Now 
we proceed to estimate in more detail the effects of the in- 
ternal magnetic field of vortices on their interactions and 
scaling. If the vortices are restricted to a two-dimensional 
film and the magnetic field can extend into the space 
above and below the film, then the interaction at large 
distances is of the form F(r) oc l/r 2 [l3j] . 

For a three-dimensional block of superconducting ma- 
terial, where the interactions occur entirely within the 
block, the magnetic field falls off over a length scale given 
by the penetration depth A. The energy of interaction 
between vortex lines in this case is 



E(r) 



8tt 2 A 2 



(2) 



where K is the zeroth-order Hankel function and <S? 
is the magnetic flux per vortex. Kq behaves logarith- 
mically as r — ► and as exp(— r/A)/r 1 / 2 as r — ► oo. 
The corresponding asymptotic forms for the force F are 
F(r) ~ (r/A) -1 as r — > and F{r) ~ exp(— r/A)/r 1 / 2 as 
r — > oo. 

We separate out the interactions into nearest-neighbor 
interactions and distant interactions. The nearest- 
neighbor interactions dominate the expansion of the vor- 
tex array with time. In the zero-field case, the distant 
interactions never become important as it is possible for 
vortices to leave the system at the walls, because the to- 
tal system pressure is zero. However, a non-zero applied 
field corresponds to a net pressure for vortices with the 
same field direction, and this constrains the maximum 
separation between vortices and leads to a stable vortex 
lattice at long time. 

For the long-range limit, the exponential decay of the 
force dominates, and we have: 



/re 1 * ± iy/netf^y/r) oc t 



(3) 



The first term of this can be approximately inverted as 
r — > oo to yield 



r ~ ln(t)(l - 



ln(lnft)) . 
2 ln(t) + 1 ' 



(4) 



At long times, this gives r oc ln(t). 

In the case of a purely two-dimensional system (2D), 
the results are the same (i 1 / 2 and ln(i) scalings). How- 
ever, in the case of a two-dimensional slab with three- 
dimensional fields (2DF), we obtain r oc t 1 ^. 

These predictions are fair heuristic approximations but 
do not match closely our observed scaling exponents. For 
example, in the 2DF case, rather than <j> = 1/3, the ob- 
served value is <j> = 0.294 ± 0.01. Similarly, the observed 
value in the 2D system is 0.414 ± 0.01, rather than the 
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predicted 0.5. A possible explanation for these results is 
that we are not yet in the long-time limit, due to long- 
lived influences of the initial conditions. During the early 
stages of the quench, the largest contribution to the evo- 
lution of the vortex distribution is not repulsion between 
vortices of the same topological charge, but rather annihi- 
lation between opposite vortices. That process continues 
into later times, although with reduced frequency. We 
will now explore the implications of such arguments. 



III. THE EFFECT OF INITIAL 
FLUCTUATIONS 

A. Zero external magnetic field 

We now consider the scaling of the inter- vortex spacing 
during the period of time in which vortex-vortex annihi- 
lation is stronger than the long-range repulsive forces. 
The following argument |2| predicts a scaling exponent 
(j) = 3/8 in two dimensions. 

In a two-dimensional sheet perpendicular to the vortex 
cores, the total topological charge contained within an 
area, radius R, of the sheet is related to the integral of 
the gradient of the phase of the order parameter around 
the boundary of that area: 



dr 



(5) 



For a mixture of positive and negative vortices, net ro- 
tation at the boundary will be due mostly to boundary 
charges within a distance r of the boundary: 



q{r)dAK I q(r)dA 

boundary « V^R (6) 



This relates the total charge fluctuations in the interior: 



which implies 



p(R) oc R 



-3/2 



(7) 



(8) 



With diffusive coarsening of domains, the domain radius 
will scale as y/t, giving p(t) oc t~ 3 / 4 . The spacing be- 
tween vortices is proportional to p^ 1 ^ 2 , leading to an 
inter- vortex scaling of i 3 / 8 in 2D. This argument applies 
equally well in a purely 2D system and a purely 3D sys- 
tem, but the growth of the domain size will be different 
in a 2DF system. 

An analysis by Jang, et al. |lfll | extends this argument 
by considering how the repulsive force between vortices 
scales with domain size, thus introducing corrections to 
the simple assumption of diffusive scaling used above. 



For the case when the inter- vortex force is proportional 
to r _1 , the force on the boundaries of the domain scale 
as 



F (x r- 1 (x yfp <x iT 3/4 

Assuming overdamped dynamics of the domain size: 
dR 



(9) 



dt 



£-3/4 



(10) 



we solve for R, with the results that R oc t 4 / 7 , p oc i~ 6 / 7 , 
and 4> — 3/7. 

A further extension of this argument was proposed by 
Liu, et al. 11] for systems of arbitrary dimension (rel- 
ative to the defect dimension), with the general result 
that 



2cj) = v 



d(d+l) 
d 2 + 2d-l 



(11) 



We can use the argument due to [10( to determine the 
vortex-annihilation scaling in the two-dimensional slab 
system, in which the domain scaling is driven by a r~ 2 
force. The result of this is that R oc i 2//5 and r oc i 3 / 10 . 

For the exponentially decaying magnetic interaction of 
the bulk 3D superconductor, there is expected to be a 
logarithmic scaling when the magnetic interaction term 
dominates the dynamics. 



B. External magnetic field 

If the external field is nonzero, there is an equilibrium 
spacing between the vortices. To treat this case, we posit 
that the interaction with the walls and remainder of the 
vortex lattice is of the same form as the magnetic inter- 
action, with a vortex placed at a distance 2R. We can 
make an argument from the vortex dynamics model to 
determine the inter-vortex spacing as a function of time. 
The methods involving initial fluctuations and vortex an- 
nihilation are more difficult to apply in this case as one 
does not simply have expanding domains, but rather the 
relaxation into a vortex lattice state with nonzero equi- 
librium topological charge. 

We assume the existence of the lattice, and look at the 
dynamics of a single vortex moving towards its equilib- 
rium position under the influence of the force 



F(r) = F (X(e 



-r/X _ -(R-r)/X 



)) 



We obtain: 



^) = -±Aln(i±44' 

K> 2 y l-a(t)- 



(12) 



(13) 



where 



a(t) = e - 2 ^-*°) e ~ R/2X , f3 = F X (14) 
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System 


Vort Dyn 


Magnetic 


Flucts Measured 


2D 


1/2 


In 


3/7 


2DF 


1/2 


1/3 


3/10 0.294 ±0.01 


3D 


1/2 


In 


3/7 0.414 ±0.01,ln 



TABLE II: Summary of predictions of vortex scaling and re- 
sults 



The equilibrium vortex spacing should scale as R = 
a/(Vh\) where a = X^/Wq. This is simply a consequence 
of distributing the flux through the sample into discrete 
vortices each carrying flux $o. We combine this with our 
predicted form, with the result that: 



(15) 



This form neglects effects due to the interaction of vor- 
tices through the order parameter, and so for very low 
fields it is expected to fail in the regime in which one 
observes power-law dynamics. This can be remedied in 
principle by adding a term with the proper force law, but 
is undesirable as it introduces more adjustable parame- 
ters. It is also expected to fail at times t < 10 during 
which the quench dynamics dominate vortex dynamics, 
but should be valid at the t — > oo limit as it reproduces 
the equilibrium vortex spacing. 

A different way of handling this problem is to proceed 
under the hypothesis that the entire effect of the external 
field is to introduce a crossover between two length scales 
- the prediction from the previous analysis for the zero 
field case and the equilibrium vortex spacing. 

The scaling in the absence of a magnetic field defines a 
length scale 1$ oc t^, whereas the Abrikosov vortex lattice 
spacing produced by the externally applied magnetic field 
is a second length scale lh oc 1/Vh relevant as t — > oo. 
The inter-vortex spacing is expected to have the form 



r(t) = ^/(i- /i~ 1/2 ) 



(16) 



such that f(x) oc x when x — > [t — ► oo, h ^ 0) and 
f(x) oc constant as x — > oo (h — ► 0, t ^ 0). 

A summary of our predictions and the results reported 
below for our numerical calculations is given in Table ITT1 



IV. SIMULATIONS 



We test this predicted form for the inter- vortex separa- 
tion by simulating the complex Ginzburg-Landau equa- 
tions coupled to a vector potential 0: 



where the covariant derivative is given by 
D = (V - iA) 



(17) 



(18) 



We must also simultaneously solve Maxwell's equations 
in the presence of a spatially varying conductivity. Using 
Maxwell's equations, we obtain |l5|: 



f) A 
at 



(19) 



We have chosen the guage V • A = to simplify Eq. El 
and Eqn. ^3 As a consequence, this must be explicitly 
satisfied by the numerical method. By using the link- 
variable method we can satisfy this constraint automati- 
cally. In this method, the vector potential is stored along 
the links between grid centers, and the order parameter 
is stored at the grid corners. Differential operators (the 
covariant derivative in particular) are evaluated as inte- 
grals around loops containing the corner of interest. The 
form of the integrals ensures that even if the gauge field 
contains a divergence, it has no effect on the numerical 
evolution. This technique has been used to simulate both 
Type-I and Type-II superconductors 0, |j| . 

Using the discretization of 0], the link variable Ui is 
defined as 



Ui 



iAiAx 



(20) 



Then, the covariant derivatives in the forward and back- 
ward direction are: 



Ui, X tP; 



x-\-i 



1px 



Aa 



j>x - Ulx-j^a 

Ax 



(21) 



(22) 



Whenever D 2 appears in the equations of motion, it is 
taken as D~D + to ensure rotational invariance. The 
covariant derivative in the superconducting current term 
is taken to be D + . This method is equivalent to making 
local gauge transformations to remove the gauge field 
from the covariant derivative, but as a result introducing 
local rotations of the phase of tp. 

Often the boundary conditions are the most compli- 
cated element of such simulations. We can simplify the 
boundary conditions by allowing there to be a thin shell 
of insulator around our sample (as we do not have peri- 
odic boundary conditions, this is consistent). The bound- 
ary conditions for the insulator are simpler as one does 
not need to consider the behavior of the superconduct- 
ing current. The latter would involve the boundary con- 
dition that W n tp — iA n 17], and because this contains 
both values of ip and A it is not very convenient to solve. 
We proceed to generate a layer of insulator by having 
a spatially-varying value of a in the Ginzburg-Landau 
equation. The value a — 1 corresponds to the supercon- 
ducting phase being stable (at zero field), whereas the 
value a = — 1 corresponds to the insulator being sta- 
ble. The order parameter will fall off to zero smoothly 
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on a scale of the order of the coherence length across the 
sharp a boundary, and so the boundary conditions are 
automatically taken into account. 

To apply an external field, we must choose the form 
of the vector potential along the boundaries. We pick a 
vector potential of the form A = H(—(y — y center )/2, (x — 
^center) /2, 0). This can pose a problem, because choosing 
a particular form of the vector potential on the bound- 
aries is equivalent to a choice of gauge. So as a conse- 
quence, we break the gauge invariance and thus break the 
translational symmetry of the system. We have not found 
an elegant way around this, but the range at which the 
effects become significant can be estimated. The break- 
ing of symmetry occurs because the gradient of the phase 
of the order parameter is slaved to the vector potential. 
The magnitude of the vector potential varies over space, 
attaining a minimum at a point which we choose when we 
apply the external field. There is a rotation of the phase 
along curves surrounding this point. If the wavelength 
of the phase rotation is smaller than the size of a grid 
cell, then artifacts will appear in the form of preferred 
locations for vortices (pinning). 

For a field H applied to a sample with cross-sectional 
area A, there is a flux cf> = HA through the sample. This 
creates <^/<E>o vortices, so the phase of the order parameter 
rotates by 2ttHA/&q along the outer boundary of the 
system, which has a length \J~A~. From the constraint for 
the physical accuracy of our system, that d x 9 < n/Ax 
where Aa; is the lattice spacing, we obtain: 



100 



Ax < 



2Hy/A 



(23) 



which constrains our lattice spacing to prevent a fail- 
ure to resolve the order parameter at the boundaries. 
This requires us to have at least N grid points for a 
two-dimensional system of edge length L, where N = 
4L 4 iJ 2 /$Q. For a three-dimensional system, this con- 
straint does not exist along the direction of the magnetic 
field, only in the plane perpendicular to it, so the cost 
there is not necessarily as bad as 0(L 6 ). Also, for suffi- 
ciently small applied magnetic flux, this constraint on the 
number of computational cells necessary is weaker than 
the constraint that we must be able to resolve the vortex 
cores: Aa: < r core . This gives the following forms for 
the number of computational cells in 2 and 3 dimensions 
(N 2D and N 3D ): 



N 2D = 



N aD = 



L' 2 



(24) 



(25) 



An 0(L 6 ) cost is simply intractable, so we must limit 
ourselves to weak applied fields. 
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FIG. 1: (Color online) This figure shows the power-law scal- 
ing of the 2DF (256 x 256 x 32 with k = 1) and 3D (92 3 with 
k = 3) systems at H — and the deviation from scaling at 
long times. The error bars were determined by the standard 
deviation between different runs. 



If the system is not in the equilibrium vortex-lattice 
state — for example, if all the vortices are in the center 
of the system, the worst-case scenario — then at a dis- 
tance r from the center the phase must rotate with a 
rate d x 6 — HA/<fi c r. There will be significant nonphysi- 
cal behavior unless r is smaller than the vortex core size 
(so that the order parameter goes to zero). Thus we 



have r r 



< r„ 



where r r 



HAAx/n4> c . This is the 



other mode of failure. In practice, this is not as strin- 
gent a constraint as there will usually be vortices spread 
evenly throughout the system during the simulation of a 
quench, which reduces the rotation of the phase near the 
center. 



We have considered ways of removing this constraint: 
for example, by evolving the gradient and magnitude 
of the order parameter, but these techniques introduce 
other technical problems. Specifically, the equations of 
motion of the gradient of the order parameter have sin- 
gularities at vortex cores, and the constraint that vortices 
be quantized must be administered separately, requiring 
that the integral of the gradient of the order parameter 
around a closed curve is 27m. 



6 




t exp(-a/h 1/2 ) 



FIG. 2: (Color online) Data collapse of the scalings at fields 
h = 0.1 to h — 0.5. The inset compares the inter- vortex 
spacing to the functional form of Eqn. 1131 The values a = 
1.626 and A = 1.810 were used to produce this collapse. These 
data are from a 3D simulation at size 92 3 and K — 3. 



V. RESULTS 



A. Three dimensional superconductor in zero 
external field 



For the three dimensional block superconductor, we 
simulated thirty systems with 96 3 grid cells with physical 
dimensions (32r core ) 3 and k = 3, starting from a random- 
ized order parameter and allowed to evolve in time. We 
also simulated five runs of a 128 3 system for longer times 
in order to observe deflections from the small-r scaling. 

The initial state of the order parameter in each cell 
was chosen by generating a random angle between and 
2tt and a random amplitude between 1 and 0, both with 
a flat distribution. In the case where the superconductor 
was restricted to a subset of the space (that is, the 2DF 
system) the amplitude was initialized at zero outside of 
the thin plate of superconductor. The three components 
of the vector potential were set to random values between 
±1 x 10~ 6 in order to break the symmetry of the field 
and reduce potential artifacts. 

We compute the average distance between vortices in- 
directly by measuring the total volume of vortices in the 
system and dividing by the cross-sectional area of a vor- 
tex, then using that density to obtain the inter-vortex 
spacing: 




(26) 



\l Jdv(i-\m\ 2 ) 

TablelTTlsummarizes our numerical results and the predic- 



10 20 30 40 50 60 70 80 90 100 
th -0.92 

FIG. 3: (Color online) Data collapse of the nonzero field scal- 
ings according to Eqn. 1161 Data are from the same simulation 
as Fig. H 



tions made here. For comparison, we provide a list |IJ of 
some previous numerical results that have been obtained 
in simulations of similar systems. The majority of such 
simulations have been done on 2D and 3D systems with 
vortex topological defects, but without a gauge field. 

For the 3D bulk superconductor, we observe that in the 
zero-field case after about t — 10, the behavior is power- 
law, with an exponent <fi = 0.414 ± 0.01, in comparison 
with the prediction from vortex-annihilation dynamics of 
(\> = 0.429. In the larger system we observe t 0A15 scaling, 
with the behavior deviating from this around t = 600 
(see Fig. [TJ, possibly due to finite size effects. 

/From our analysis, we also expect a region of logarith- 
mic scaling. At long times r varies logarithmically over a 
time interval determined by the system size. For all lat- 
tice sizes, the behavior between t = 5 and t = 20 seems 
to follow a logarithmic curve which then transitions to 
the power-law behavior. 



B. Two-dimensional superconductor 

For the case of a two dimensional sheet with three di- 
mensional fields, we simulated five runs of a system of 
size 256 x 256 x 32 grid cells, with k = 1. The full pre- 
diction for that case involves a region of t 3 ^ 7 scaling, fol- 
lowed by a region of £ 3 / 10 scaling for vortex-annihilation 
dominated dynamics. As a consequence, for intermedi- 
ate values of n there is no large scaling range for either 
regime without waiting long times in the simulation. The 
smaller the value of k, the shorter the time necessary to 
observe the asymptotic regime and the effect of external 
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FIG. 4: (Color online) Data collapse of the nonzero field time 
derivatives. (5 — —1.5. Data are from the same simulation as 
Fig.0 



three-dimensional fields. 

We observe a scaling of i°- 294±0 - 01 ! which is in agree- 
ment with the prediction for vortex-annihilation domi- 
nated scaling in the large-r limit. The results arc shown 
in Fig. □ 

Eventually, the separation between vortices becomes 
on the order of the system size and finite size effects ap- 
pear in the results which puts the final limit on the size 
of the scaling regime we may observe. This manifests 
itself as an increase in the fluctuations across different 
runs and as a sudden sharp upturn in r compared to 
the power-law prediction, as shown by the data for times 
t > 600 in Fig. □ 



C. Non-zero field data collapse 

The results for a non-zero applied field are in Fig. [21 
Here we also plot the best-fit from the predicted form of 
Eqn. El a l on g with each curve. The fits are satisfactory, 
but break down at short times during the initial quench. 

When we plot the data in terms of the combined vari- 
ables given by Eqn. 1151 we can observe that the time and 
field-strength variables are not independent. The family 
of curves we observe will collapse onto a single universal 
curve given the proper choice of scaling of the variables 
as we have described in section IIII Bl The behavior in 
the absence of a magnetic field also has a crossover - be- 
tween power-law scaling and logarithmic scaling. The 
model which generates this collapse does not take into 
account the power-law scaling and so we expect devia- 
tions from the collapse in the case that there is a large 
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FIG. 5: (Color online) Error landscape for the r data collapse. 



regime of power-law scaling Qi — > 0). Deviations from 
the collapse in our simulation are observed at times less 
than t = 20 due to the initial quench behavior which is 
not captured by this model. At sufficiently small fields 
that the Abrikosov lattice spacing is comparable to the 
system size we also expect a failure of the collapse at 
long times - in the zero field case for the 92 3 system this 
occurs at t > 600. 

In order to capture the power-law scaling behavior, we 
propose a different form of the data collapse given in Eqn. 
UHl Plotting the data accordingly we obtain Fig. The 
scaling exponent <fi was adjusted to give the best visual 
data collapse, which occurred for <j> = 0.46 ± 0.2. 

The behavior at short times does not collapse, which is 
consistent with the concept behind this scaling: at short 
times, the scaling is dominated neither by vortex inter- 
actions nor the external magnetic field. Instead, it is the 
local dynamics of the order parameter which contribute 
to the inter-vortex spacing. During this time, the mea- 
sure we have used to determine the inter-vortex spacing 
cannot be considered meaningful as there are not yet dis- 
tinct vortex cores. 

We have also examined the time derivative of the vor- 
tex spacing, as this is the basis of the force-law arguments 
to determine vortex scaling. The behavior of f as a func- 
tion of r at nonzero external fields asymptotes to zero at 
values of r determined by the external field, i.e. the equi- 
librium vortex spacings. Rescaling r to map these to the 
same point results in a new variable r' — \fhr. However, 
if r scales as r 13 in the limit of small r, then we must also 
rescale r as follows: r' = h^' 2 r. This collapses the data 
for f onto a single curve, as seen in Fig. 0] The value of 
(3 which best collapses the curves is sa —1.5, which corre- 
sponds to a value of the scaling exponent <f> = 0.4. Unlike 
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FIG. 6: (Color online) Equilibrium vortex spacing at nonzero 
field. The line is a linear guide. Data are from the same 
simulation as Fig. |2] 
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h l/24> 

FIG. 7: (Color online) Scaling of the vortex spacing in the 
limit of short times. <f> = 0.46. The line is a linear guide. 
Data are from the same simulation as Fig. |5] 



the r data collapse, the error landscape of the f collapse 
has a well-defined minimum at finite <f>, although there is 
still a gradual decrease in the error as 4> becomes unphys- 
ically large. Figure [S] shows the error in the data collapse 
measured by A = V_\. fc ~ Vk{ x i)) 2 where ir, is the 

ith rescaled independent variable data point and yj and 
yk are the interpolated values of the rescaled dependant 
variable for the jth and feth field strengths. 

The data collapse of Fig. is actually an autonomous 
counterpart to that of Eqn. 1161 as can be readily verified 
by elimination of t. 

Now we turn to the small field data, which seems not to 
obey the data collapse so cleanly. For the r scaling they 
do not fall onto the data collapse lines nearly as well as 
the higher fields. However, for the f scaling this failure of 
data collapse is not nearly as pronounced. This may be 
associated with a bias in the initial spacing introduced 
by the period immediately after the quench. That is, 
for the r scaling it may be that rather than measuring 
time from the simulation start, we must measure it from 
the end of the quench period when vortex interactions 
become the dominant effect. However, for the r scaling 
we avoid this by eliminating the time variable and thus 
removing the problem of an arbitrary initial offset. We 
will now examine other predictions of the data collapse 
in order to understand the small field behavior. 

The length scale we observe at long times corresponds 
to the Abrikosov vortex lattice spacing. Finite system 
size effects can potentially produce a crossover to the 
system scale and are more likely to do so for weak field 
values. The value of the spacing at the longest times 



accessed by our simulations is plotted versus the field 
strength in Fig. For sufficiently large fields h > 0.2 
(which corresponds to hr 1 / 2 < 2.24) we see a conver- 
gence to the predicted scaling within the time range of 
our simulation. For small fields we never quite reach the 
equilibrium scaling. We have however also looked at a 
128 3 system at low fields and during the time range of 
our simulations (up to t = 400 in this case) we did not 
observe a significant departure from the 92 3 curves, so 
the slight deviation from Abrikosov scaling is more likely 
to be a finite time effect than a finite size effect. 

We can also make a prediction as to the scaling of 
r(h) at a fixed small time t\. Inserting into the scaling 
form, we get r(h,h) = fty-T^hr 1 / 2 ). We wish to 
Taylor expand / around t = 0, which we may safely do 
after writing / = with the result that r(h,ti) — 

tr^Khh 1 / 2 ^) aCi + c a H 1 / 2 +t\-* + o{h 2 -^). 

We know that C\ = due to the required asymptotic 
behavior of the scaling function, meaning that we can 
expect r{H, t\) oc /i 1 / 2 ^ for fixed small times. The agree- 
ment with this prediction is satisfactory, as shown in Fig. 
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VI. CONCLUSIONS 

We observed both the case of a three dimensional block 
of superconductor and a two dimensional sheet via nu- 
merical simulation. In the three dimensional block su- 
perconductor, we also examined the effect of nonzero ex- 
ternal fields, and predicted a scaling form based on the 
forces between vortices. With this scaling form, we ob- 
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tained data collapse to a universal curve independent of 
the value of the external magnetic field. At zero exter- 
nal field, the dynamical critical exponent 4> was found to 
be in agreement with fluctuation/annihilation arguments 
rather than the pure vortex-dynamical argument in both 
the 3D and 2DF systems - that is, in the 3D case we mea- 
sured (f> = 0.414 ± 0.01 in comparison with the predicted 
3/7 and in the 2DF case we measured <j> = 0.294 ± 0.01 
in comparison with the predicted 3/10. We also ob- 
served data suggestive of the logarithmic scaling regime 
at short times between the recovery of the order parame- 
ter from the quench and the power-law scaling regime. At 
long times, we observed deviation from power-law scal- 
ing which appears to be a finite size effect. We never 
observed a regime in which all vortices remaining were 
of one direction/magnitude of topological charge in the 
zero-field case. However, with a nonzero field, the major- 
ity of vortices in the system are aligned with that field, 
and so the argument based on vortex dynamics seemed 
to be more able to predict the scaling form than in the 
zero-field cases, in which vortex annihilation was always a 
significant contribution. We have proposed a form for the 
behavior of the vortex spacing in the presence of external 
magnetic fields which collapses the data onto a universal 
curve. We found better collapse when this was applied 
to the force law as opposed to the actual vortex spacing 
likely due to the initial time offset in which the order 
parameter recovers from the quench. This allowed us to 
explain the behavior at nonzero magnetic field in terms 
of a crossover from the zero-field scaling to a fixed length 
scale given by the Abrikosov vortex lattice spacing. 

We would like to thank Arttu Rajantie, Badri Athreya 
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